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Abstract. We study the order-disorder transition in colloidal suspensions under 
shear flow by performing Brownian dynamics simulations. We characterize the 
transition in terms of a statistical property of time-dependent maximum value of the 
structure factor. We find that its power spectrum exhibits the power-law behaviour 
only in the ordered phase. The power-law exponent is approximately -2 at frequencies 
greater than the magnitude of the shear rate, while the power spectrum exhibits the 
l//-type fluctuations in the lower frequency regime. 
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1. Introduction 

Order-disorder transitions such as solid-liquid transitions and ferromagnetic transitions 
are distinctive phenomena in equilibrium systems. The understanding of the transitions 
has been a cornerstone for developments in equilibrium statistical mechanics. Here, 
the notion of "order" is not restricted to equilibrium systems. For example, one can 
modify a system exhibiting the order-disorder transition under equilibrium conditions 
so that it can be observed even under non-equilibrium conditions. Since non- 
equilibrium statistical mechanics has not been established as yet, there is no systematic 
understanding of this type of transition from the viewpoint of statistical mechanics. 
Therefore, it is important to study a typical example related to this question. 

As a simple and realistic example, we consider a system consisting of colloidal 
particles suspended in a liquid. It is known that this system exhibits an order- disorder 
transition under equilibrium conditions. In particular, the so-called "colloidal crystal" is 
observed in the ordered phase. This system may be regarded as an ideal model system 
because one can observe such a crystalline structure by using a microscope and also 
because one can control the system as desired. 

Recently, a couple of experiments have been reported on colloidal suspensions 
under shear flow [TJ [2]. Holmqvist et al. obtained the phase diagram for colloidal 
suspensions under stationary Couette flow. They observed Bragg reflections of laser 
light and measured the Bragg peak intensity of the first Debye-Scherrer ring. From the 
time dependence of the total Bragg peak intensity, they determined the crystal growth 
rate and induction time. The phase boundary in their phase diagram is determined by 
extrapolating the growth rate and induction time are extrapolated to zero. 

Furthermore, there are several reports on numerical experiments of colloidal 
suspensions under shear flow. Butler and Harrowell performed Brownian dynamics 
simulations of colloidal particles under shear flow (3JH]. They obtained a phase diagram 
of temperature versus shear rate. The phase boundary of their phase diagram was 
determined from the long-time average of the intensity of the Bragg peak with the 
wavevector aligned in the shear gradient direction. The 'order' in their paper was 
assumed to appear when the intensity exceeded half of the intensity value for scattering 
from a body-centred cubic crystal aligned in the direction of the velocity gradient. 

These results suggest the existence of the order- disorder transition. Here, let us 
recall that crystal is defined as the state with a translational symmetry breaking. Thus, 
if we attempt to determine the 'crystal' phase, we should investigate the structure factor 
S(k) for systems under shear flow. However, from a simple consideration, we find that 
there is no crystal phase in the rigorous sense, even if the existence of an ordered state 
is suggested by several measurements. Now, the objective of this study is to find a 
useful and thorough characterization of the ordered state for systems under shear flow. 
We attempt this characterization by performing Brownian dynamics simulations of the 
colloidal suspensions under shear flow. 

This paper is organized as follows. In section 2, the model that is investigated 
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is introduced. Further, before discussing non-equilibrium cases in section 3 as 
preliminaries, we review the order- disorder (crystal-fluid) transition under equilibrium 
conditions and confirm the conventional criteria for the crystallization. Section 3 
comprises the main part of this paper. We characterize the ordered phase in terms 
of the time-dependent quantity s m (t), which is defined as the first maximum of the 
structure factor s(k,t) for configuration of the particles at time t. We find that S m (u), 
which is the power spectrum of s m (t), exhibits a clear transition when the temperature 
is changed. Concretely, in the ordered phase, the power spectrum exhibits the power- 
law behaviour, and there are two different power-law exponent regimes divided by the 
crossover frequency that is determined by the shear rate 7: S m (uj) ~ uj~ 2 in the high 
frequency regime (a; > 2717) and l//-type fluctuations in the low frequency regime 
(ui < 271-7). Meanwhile, S m (uj) in the disordered phase is similar to the white noise 
spectrum. Future problems are presented in the final section. 



2. Preliminaries 



We consider a system consisting of N colloidal particles suspended in a fluid where the 
stationary planar shear flow is realized. The system is confined to a cubic cell with a 
length L. The x-axis and z-axis are chosen to be the directions of the shear velocity and 
velocity gradient, respectively. We impose Lees-Edwards periodic boundary conditions 
[5, 6J to avoid peculiarities near the boundaries of the cell. 

We assume that Langevin dynamics can describe the motion of the colloidal 
particles. Concretely, the force exerted from the fluid is represented by the Stokes force 
and Gaussian noise. In other words, we neglect the so-called hydrodynamic effects. 
Then, the particle positions 7*j(t) = (xi(t), yi(t), Zi(t)), where 1 < i < N, obey the 
Langevin equations 

l^t = ~ E W (l r * - r iD + + C 1 ) 

where rj is a friction coefficient; 7, the shear rate; and e x , the unit vector that is parallel 
to the a:-axis. The variable = {^f(t),^(t),^(t)) represents the Gaussian noise that 
satisfies 

[m)^(t')) = 2 v k B T6 tJ 6^6(t - t'). (2) 

Here, &b is the Boltzmann constant and T is the temperature. The superscripts a and 
(3 represent the Cartesian components. Each pair of particles interacts via a screened 
Yukawa potential 

f rr ( exp(-^(r - a)) exp(-/t(r c - a)) \ 

[ 0, otherwise. 

Here, r is the distance between the particles; r c , the cutoff length that simplifies the 
calculation in numerical simulations; and k, a Debye screening parameter. 
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In this study, all the quantities are converted to dimensionless forms by setting 
o = U = r//k B T = k B = 1. We estimate the correspondence between these 
parameters and those of real experimental systems to be a ~ 10 2 nm, Uq ~ lOk-^T 
and k-QT/r] = 10 _11 m 2 s _1 . In our simulation, we assume that kg = 5.8, L/a = 10, 
r c = 2.5 and Na 3 /L 3 = 1. In other words, the Debye screening length k~ x corresponds 
to 17 nm. The typical values of the parameters used in our simulations are T ~ 0.1 
and 7 = 0.001, and this situation corresponds to experimental systems wherein the 
temperature is 300 K and the shear rate is 1 s -1 . The systems that have the values are 
available by laboratory experiments. 

In our simulations, we discretized (CD) with the time step At = 0.0025. Note that 
in the arguments below, (■ ■ •) represents the statistical average in the steady states. In 
the calculation we performed, we estimated (A) to be J^ +T dtA(t)/r, where to and r 
were chosen as that larger than 10 3 , because we had confirmed that the relaxation time 
is approximately 10 2 . 

Before considering the behaviour of colloidal suspensions under shear flow, we 
review the transition observed in the system under the equilibrium condition (7 = 0). 
It is expected that the crystalline arrangement of colloidal particles is observed in the 
ordered phase of this system, while they acquire a random configuration in the disordered 
phase. In order to detect the order-disorder transition, one relies on the definition of a 
crystal according to which the statistical weight for configurations of colloidal particles 
breaks translational and rotational symmetries. The translational symmetry breaking 
can be quantified by S(k), which is defined as 

s(k) = (^- r d<p r ae sine hp{k)\ 2 \ . (4) 



Here, the angles e and <ft are defined as k = (/ccos</>sin#, /csin0sin#, fccos#) with 
k = \k\, and p(k) is the Fourier transform of the number density p(r): 



00 

3 



p(k) = / d r p(r) exp(— ik ■ r) (5) 



—00 



with 

1 - 

^) = ^E<^-^)- ( 6 ) 

i=l 

When S(k) has a component that is expressed by Dirac's delta function, the system 
is assumed to exhibit the translational symmetry breaking. Note that S(k) can be 
measured experimentally because it is related to the Bragg peak intensity of the laser 
light scattering. 

In our simulations, instead of S(k), we use Si(k) defined by 

n /in •< a f L ^ 2 1 // / r im a \ sin (At) sin(27rr/L) 

S L (k) = 1 + 4vrp dr ((g(r; {r,})) - 1) ^ ^ ^ ; , (7) 
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where g(r; {ri}) represents the radial distribution function for a given configuration 
{ri}, that is, 

1 dn(r) 



g(r; {r t }) 



(8) 



Arrpr 2 dr 

Here, dn(r) is the average number of particles at distances between r and r + dr from 
any particle, and its average is taken over all the particles. Note that (g(r; {ri})) is 
equal to the standard radial distribution function g{r). It should be noted that g{r) is 
defined in the range where r < L/2, owing to the periodic boundary conditions. The 
term sin(27rr/L)/(27rr/L) appended in the integrand is the window function to reduce 
the termination effects resulting from the finite upper limit 018]. Note that S^k), as 
defined above, approaches S(k) in the thermodynamic limit L — > oo. 





Figure 1. Plot of g{r) against r/a. 
T = 0.12 (O : open circles) and 
0.18 (*: stars) under equilibrium 
conditions. 



Figure 2. Plot of SL{k) against 
kcr. T = 0.12 (O : open 
circles) and 0.18 (*: stars) under 
equilibrium conditions.. 



As observed in figure [2l we cannot observe infinitely sharp peaks owing to finite 
size effects. However, it is known that there are empirical criteria for detecting the 
transition; one of these criteria is Hansen-Verlet's rule [9j [TO]. According to this rule, 
a fluid freezes when the first maximum value of the static structure factor S m exceeds 
2.85. This criterion has been tested and validated for various systems [SHU]- Indeed, in 
our system, we find that S m exhibits a discontinuous jump at the transition temperature 
T c whose value is estimated between 0.16 and 0.165, as shown in figureEJ Further, when 
the temperature T is lower than the transition temperature T c , S m exceeds 2.85, while 
S m is less than 2.85 in the higher temperature regime. Thus, we conclude that the 
order-disorder transition is observed. 

As another order parameter for indicating the phase transition, we consider the 
bond-orientational order parameter Q 6 [121 [13] , which is determined by the set of bond 
vectors {fj} as follows: 

1 /2 

«« = ((|ti«i 2 ) )■ w 
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with 



1 N h 



(10) 



Here, each bond vector corresponds to the relative vector between the neighbouring 
particles; iVb is the number of bond vectors; Yq, the spherical harmonics function of 
degree six; and 0(ri) and <p{ri), the polar and azimuthal angles of fj, respectively. 
Here, we have defined the neighbouring particles for a given particle as those within the 
sphere of radius r c around the given particle, where r c is chosen as the first minimum 
of g{r). The quantity Qq represents the degree of breakage of the continuous rotational 
symmetries, particularly, the 6-fold rotational symmetry of the configuration of the 
particles. Its value is 0.57452 for a face-centred cubic (fee) crystal, 0.51069 for a body- 
centred cubic (bec) crystal and for liquids. We show the temperature dependence of 
Qq in figure HI It is observed that Qq decreases between T = 0.16 and T = 0.165 in a 
discontinuous manner. This result is consistent with that indicated by Hansen-Verlet's 
rule. 
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Figure 3. S m as a function of 
T under the equilibrium condition 
7 = 0. The average values with 
error bars are displayed using 10 
samples for each temperature. 



Figure 4. Qq as a function of 
T under the equilibrium condition 
7 = 0. The average values with 
error bars are displayed using 10 
samples for each temperature. 



3. Question and Result 

Even for colloidal suspensions under shear flow, we can measure S(k) and Qq in a manner 
identical to that of the equilibrium cases. The results with 7 = 0.001 are displayed in 
figures [5] and [6j These graphs are similar to those for the equilibrium system. In the 
temperature regime lower than T = 0.16, S m exceeds 2.85, while it dips from 2.85 in 
the regime higher than T = 0.17. Apparently, this result indicates the existence of an 
order-disorder transition in this situation. Similarly, Qq shows a clear difference between 
the low temperature regime and the high temperature regime, as shown in figure [6j 
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Figure 5. S m as a function of T at 
7 = 0.001. The average values with 
the error bars are displayed using 
10 samples for each temperature. 



Figure 6. Qg as a function of T at 
7 = 0.001. The average values with 
the error bars are displayed using 
10 samples for each temperature. 



We recall that a crystal is defined as the state that breaks the translational 
symmetry of the statistical weight for the configurations of the particles. In order 
to simplify the above argument, we first consider the case T = where particles move 
coherently and form two-dimensional crystals in (x, y) planes for a fixed z. In this state, 
the translational symmetries in the x and y directions are broken for a given z. However, 
the spatial period in the z direction is time dependent. Thus, by considering the 
ensembles generated by the time evolution, we expect that the translational symmetry 
in the z direction recovers. Next, we consider the finite temperature cases. A 
translational symmetry is expected to exist in the z direction, and according to Mermin- 
Wagner's theorem, which states that there is no translational symmetry breaking in two- 
dimensional systems, we do not expect the symmetry breaking to occur in the (x, y) 
planes. Thus, we conclude that there is no crystal in colloidal suspensions under shear 
flow. 

Although the structure factor does not involve Dirac's delta peak even in the 
thermodynamics limit, figures [5] and [6] suggest the existence of a discontinuous transition. 
Therefore, the nature of this transition should be different from that in equilibrium 
systems. In order to further clarify the above difference, we attempt to study this 
transition from another viewpoint. 

In particular, we consider the dynamical features of non-equilibrium systems 
possess. First, let us observe the variation of the structure factor with time. As an 
example, we define the time-dependent structure factor s(k,t) that is determined for 
each particle configuration at time t: 

s . , f R , / , r , , \ sin (At) sin(27rr/L) , . 

s(k,t) = l + 4irpj dr(g(r;{n(t)})-l) ^ ' ^ ' 1 . (11) 

Note that S^k) defined in (j7j) is equal to (s(k, t)). In figure [TJ the results for two 
different cases T = 0.14 and 0.18 are displayed by fixing the shear rate to be 7 = 0.001. 
From this figure, it is observed that the maximum intensity of s(k,t) at T = 0.14 
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varies with time more significantly than that at T = 0.18. Since there appears to be a 
qualitative difference, we focus on the time dependence of the first maximum of s(k,t) 
with respect to k, which is denoted by s m (t). 
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Figure 7. Time evolution of the structure factor for T — 0.14 (left) and T = 0.18 
(right). 




Figure 8. s m (t) for T = 0.14 (O : open circles) and T = 0.18 (*: stars). 



In figure [HI typical data of s m (t) are displayed. It is clearly observed that s m (t) 
exhibits a considerably larger fluctuation in the low temperature case (T = 0.14) as 
compared to that in the high temperature case (T = 0.18). In order to characterize 
the difference between the two cases quantitatively, we consider the spectra of the 
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fluctuations, which are defined by 



dt s m (t) exp(—iut) 



(12) 



The shapes of the power spectra shown in figure [9] indicate a distinct transition at a 
certain temperature T c . Indeed, for the temperatures T = 0.14 and 0.16, the power- 



law behaviour S„ 
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is observed in the frequency regime u > 27ry. Moreover, 



focusing on the behaviour in the low frequency regime u < 2tt a /, we observe the 1//- 
type fluctuation. Meanwhile, the spectrum becomes flat at the temperature T = 0.18. 
Therefore, for example, by plotting S* m (271-7) against T, we can observe a discontinuous 
transition at T = T c , as shown in figure [TUl 




co/2jiy 



Figure 9. Spectra S m (u>) as a function of = ui/2tt a {. T = 0.14 (closed diamonds), 
0.16 (open circles) and 0.18 (stars). The solid line represents the w -2 slope, and the 
dashed line denotes the w _1 slope. 
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Figure 10. SVn (271-7) as a function of T. The average values with error bars are 
determined by using 10 samples for each temperature. 



4. Concluding Remarks 

The results in this paper motivate us to study the system in more detail. Before 
concluding the paper, we address two important future problems. First, the mechanism 
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of the power-law behaviours of S m (uj) should be elucidated on the basis of "defects" in 
the ordered phase. After defining the defects in a suitable manner, we may quantitatively 
characterize the destruction of the crystal-like structure and its subsequent restoration 
with time. It is natural to conjecture that the power- law behaviour S m (uj) ~ lo~ 2 might 
be related to the generation and annihilation of defects in the ordered phase under 
shear flow. Therefore, by focusing on the elementary processes of defect dynamics, 
we may understand this power-law behaviour. Furthermore, the power-law behaviour 
S m (uj) ~ uo~ l observed in the low frequency regime suggests the existence of a more 
complex mechanism. As one possibility, a cooperative phenomenon involving the 
'defects' may occur. In order to explore this possibility, we should investigate the spatial 
correlation of the dynamical events of the defects. 

The second problem is related to the thermodynamic aspects of the order-disorder 
transition in colloidal dispersions under shear flow. From the analogy of the order- 
disorder transition in equilibrium systems, we conjecture the existence of a latent heat 
associated with the order-disorder transitions in non-equilibrium systems. However, 
the heat is apparently generated because the system is in a non-equilibrium steady 
state. With regard to this problem, Oono and Paniconi proposed a remarkable concept 
in which the excess heat plays a prominent role in a thermodynamic framework for 
non-equilibrium steady states [i4j. Indeed, by employing this quantity, the second law 
of thermodynamics has been extended to transitions between non-equilibrium steady 
states within the Langevin description [15]. Thus, we can consider the latent (excess) 
heat even for colloidal suspensions under shear flow. 

In summary, we have characterized the order- disorder transition of colloidal 
suspensions under shear flow by employing dynamical features of the structure factor 
because there exists no 'crystal' under the shear flow. Defining a new order parameter 
S(uj), we have detected the transition by measuring the temperature dependence of 
S(2Trj), as shown in figure fTUl We also have found that the discontinuous transition 
accompanies the l//-type fluctuations in the low frequency regime. By a more detailed 
study of the above-mentioned problems, we hope to gain a thorough understanding of 
the nature of non-equilibrium phase transitions. 
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